########################################
# Who wants an independent court?
# Leiras, Giraudy, Tuñón & Wagner.
# Replication file. Main program.
########################################

# set WD
setwd("~/Dropbox/TSJ provinciales/JOP/WWIC Replication files")
# loading data. 
load("JOP_WWIC_repdata_final.RData")
# loading functions
source("JOP_WWIC_functions.R")
# load libraries
require(survival)
require(foreign)
require(nnet)
#require(reshape2)
library(xtable)


############################
# Descriptive stats
############################

# Calculate mean and SD of TENURE
mean(figs.data$total_tenure)
sd(figs.data$total_tenure)

#### calculate amount of observations per type of exit
figs.data$motive.quidPol <- figs.data$motive1
table(figs.data$motive.quidPol) # non NA observations, 0 is continues in office
sum(is.na(figs.data$motive.quidPol)) # NA obs 
length(figs.data$motive.quidPol) # total observations
sum(figs.data$motive.quidPol==2 | figs.data$motive.quidPol==1 | 
      figs.data$motive.quidPol==4 , na.rm=T) # TOTAL EXITS - nonzero observations.


############################
# Descriptive figures
############################
head(polfigs.data)

# creating one dataset for pol exits
table(polfigs.data$motive1_head)
pol.exit <- na.omit(polfigs.data[polfigs.data$motive1_head==2,])
# and respective survival object
survival.long <- Surv(pol.exit$total_tenure_head)
motives.fit.long <- survfit(survival.long ~ 1)



# Figure 1. Kaplan Meier by exit. 
par(mfrow=c(1,1), oma=c(0,0, 0,0))
# generating survival object. 
survival <- with(figs.data, Surv(total_tenure, outin2008 == 1))
global.fit <- survfit(survival~1)

# general plot, all exits.
plot(global.fit, main="Judicial Survival (all exits) \n Kaplan-Meier estimates", 
     mark.time=F, lwd=3,
     conf.int = F, xlab="Judge's years in office", ylab="survival function", col="Navy Blue")

motives.fit.quidPol <- survfit(survival ~ figs.data$motive.quidPol)
plot(motives.fit.quidPol, main="Judicial Survival by Risk \n Kaplan-Meier estimates", mark.time=F, conf.int = FALSE, 
     xlab="Judge's years in office", ylab="survival function", col=c("white","Navy Blue", "Navy Blue", "Navy Blue"),
     lty=c(0,4,1,3), lwd=3)
leg.text <- c("", "Natural", 
              "Political",
              "Federal Intervention")
legend(x=5,y=1.05, bty="n", text.col=c("white","Navy Blue", "Navy Blue", "Navy Blue"), 
       legend=leg.text, cex=0.9,lty=c(0,4,1,3), lwd=3)



# Figure 2. Electoral timing of judicial purges (political exits)
plot(motives.fit.long, main="Electoral Timing of Judicial Purges (Political Exits) \n Kaplan-Meier Estimates", 
     mark.time=F, conf.int = F, xlab="Governor's years in office.", ylab="survival function", col="Navy Blue", lwd=1.75)
#lines(motives.fit.long, mark.time=F, conf.int = F, col="Navy Blue", lwd=3.5)
# lines(surv.fit.natural, mark.time=F, conf.int = F, col="red", lwd=3.5)
abline(v=4, col="gray", lty=1, lwd=3)
abline(v=8, col="gray", lty=1, lwd=3)
abline(v=12, col="gray", lty=1, lwd=3)
lines(motives.fit.long, mark.time=F, conf.int = F, col="Navy Blue", lwd=3)




# Figure 3A. Survival by province (all motives)
province.name <- c("Buenos Aires", "Ciudad de Bs As", "Catamarca", "Chaco", "Chubut", "Córdoba", 
                   "Corrientes", "Entre Rios", "Formosa", "Jujuy", "La Pampa", "La Rioja", "Mendoza", 
                   "Misiones", "Neuquen",  "Rio Negro", "Salta", "San Juan", "San Luis", "Santa Cruz", 
                   "Santa Fe", "Stgo. del Estero", "T. del Fuego", "Tucuman")

par(mfrow=c(4,6), oma=c(3,3, 4,0))
for (i in 1:24){
  province(i)
}
mtext(text="Kaplan-Meier estimate, by province", 3, cex=1.2, outer=T)
mtext(text="survival function", 2, cex=1, outer=T)
mtext(text="years", 1, cex=1, outer=T)


# Figure 3B. Survival by province (political exit)
pol.data <- figs.data[figs.data$motive.quidPol==2,]
survival2 <- with(pol.data, Surv(total_tenure, outin2008 == 1))
province.name <- c("Buenos Aires", "Ciudad de Bs As", "Catamarca", "Chaco", "Chubut", "Córdoba", 
                   "Corrientes", "Entre Rios", "Formosa", "Jujuy", "La Pampa", "La Rioja", "Mendoza", 
                   "Misiones", "Neuquen",  "Rio Negro", "Salta", "San Juan", "San Luis", "Santa Cruz", 
                   "Santa Fe", "Stgo. del Estero", "T. del Fuego", "Tucuman")

par(mfrow=c(4,6), oma=c(3,3, 4,0))
for (i in 1:24){
  province.pol(i)
}
mtext(text="Kaplan-Meier estimate, by province", 3, cex=1.2, outer=T)
mtext(text="survival function", 2, cex=1, outer=T)
mtext(text="years", 1, cex=1, outer=T)


##################################################################
##### Statistical Analysis.
##################################################################

# Before running the model, choose the level of the outcome that will be the baseline
# Here, baseline=staying in office
data$motive <- relevel(as.factor(data$motive.quidPol), ref = "0")

#### Creating necessary terms for analysis

data$fragmentation <- 1 - data$majoritydesign
data$frag.turnover <- data$fragmentation * data$turnovert2

# Polynomial terms for parametric model
data$t2 <- data$tenure^2
data$t3 <- data$tenure^3
data$t4 <- data$tenure^4


# PARAMETRIC MODEL
library(nnet)
param.main <- multinom(motive ~ tenure + t2 + t3 + fragmentation + turnovert2 
                 + frag.turnover, data = data)
param.main.complete <- multinom(motive ~ tenure + t2 + t3 + fragmentation + turnovert2 
                                + frag.turnover + term_limit + poplog + head + lemas, data = data)

param.main
param.main.complete


summary(param.main)$coefficients[2,]
summary(param.main)$standard.errors[2,]

summary(param.main.complete)$coefficients[2,]
summary(param.main.complete)$standard.errors[2,]


#function to obtain one set of simulated betas
sim.betas  <- function(coefs, se){
  sim.betas <- matrix(NA,dim(coefs)[1], dim(coefs)[2])
  for (j in 1:dim(coefs)[1]){
    for(i in 1:dim(coefs)[2]){
      sim.betas[j,i] <- rnorm(1, coefs[j,i], se[j,i])}
  }
  return(sim.betas)
}

sim.probs <- function(covariates, sims=10000, coefs, se){
  probs <- matrix(NA, sims, 8)
  for (i in seq(1, 8,1)){
    sim.data <- c(1, i, i^2, i^3, covariates)
    for (k in 1:sims){
      betas <- sim.betas(coefs=coefs, se=se)
      den <- 1 + exp(sum(betas[1,]*sim.data)) + 
        exp(sum(betas[2,]*sim.data)) + exp(sum(betas[3,]*sim.data))
      probs[k,i] <- (exp(sum(betas[2,]*sim.data)))/den
      if (round(k/100)==k/100) print(k)
    }}
  return(probs)
}

# Simulated prob of political exit
# Main model
baseline.probs <- sim.probs(covariates=c(0,0,0), 
                            coefs=summary(param.main)$coefficients, 
                            se=summary(param.main)$standard.errors) # baseline

only.fragmentation <- sim.probs(covariates=c(1,0,0), 
                           coefs=summary(param.main)$coefficients, 
                           se=summary(param.main)$standard.errors) # majority

only.turnover <- sim.probs(covariates=c(0,1,0), 
                           coefs=summary(param.main)$coefficients, 
                           se=summary(param.main)$standard.errors) # expected turnover

interaction <- sim.probs(covariates=c(1,1,1), 
                         coefs=summary(param.main)$coefficients, 
                         se=summary(param.main)$standard.errors) # majority exp turn and interaction

# getting summary stats from simulations 
plot.data <- function(matrix){
  summ <- matrix(NA, 8,3)
  for (i in 1:8){
    temp <- matrix[,i]
    summ[i,] <- c(mean(temp), quantile(temp, probs=c(.05, .95)))
  }
  colnames(summ) <- c("mean", "CI low", "CI high")
  return(summ)
}

basic.baseline <- plot.data(baseline.probs)
basic.fragmentation <- plot.data(only.fragmentation)
basic.turnover <- plot.data(only.turnover)
basic.int <- plot.data(interaction)


# Plot for paper
par(mfrow=c(1,4), oma=c(2,0,2,0), mar=c(5,5,5,2), xpd=F, cex=.85)
plot(1:8, basic.fragmentation[,1], type="l",
     ylim=c(0,.65), xlim=c(1,8), col="navyblue",lwd=3, 
     xlab="Tenure (years)", ylab="Prob. of political exit",
     main="PANEL I \n Simulated Effect of Fragmentation \n when there is no Expectation of Turnover")
lines(1:8, basic.fragmentation[,2], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.fragmentation[,3], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.baseline[,1], col="navyblue",lwd=3, lty=2)
lines(1:8, basic.baseline[,2], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.baseline[,3], col="navyblue",lwd=.7, lty=2)
legend("bottomleft", c("Fragmentation", "No Fragmentation"), 
       lty=c(1,2), lwd=2, col=c("navyblue"),
       horiz=F, bty="n", cex=.9, y.intersp=.8)

plot(1:8, basic.int[,1], type="l",
     ylim=c(0,.65), xlim=c(1,8), col="navyblue",lwd=3, 
     xlab="Tenure (years)", ylab="Prob. of political exit",
     main="PANEL II \n Simulated Effect of Fragmentation \n when there is Expectation of Turnover")
lines(1:8, basic.int[,2], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.int[,3], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.turnover[,1], col="navyblue",lwd=3, lty=2)
lines(1:8, basic.turnover[,2], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.turnover[,3], col="navyblue",lwd=.7, lty=2)
legend("bottomleft", c("Fragmentation", "No Fragmentation"), 
       lty=c(1,2), lwd=2, col=c("navyblue"),
       horiz=F, bty="n", cex=.9, y.intersp=.8)

plot(1:8, basic.baseline[,1], type="l",
     ylim=c(0,.65), xlim=c(1,8), col="navyblue",lwd=3,, lty=2, 
     xlab="Tenure (years)", ylab="Prob. of political exit",
     main="PANEL III \n Simulated Effect of Expected Turnover \n when there is no Fragmentation")
lines(1:8, basic.baseline[,2], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.baseline[,3], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.turnover[,1], col="navyblue",lwd=3)
lines(1:8, basic.turnover[,2], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.turnover[,3], col="navyblue",lwd=.7, lty=2)
legend("bottomleft", c("Expectation of Turnover","No Expectation of Turnover"), 
       lty=c(1,2), lwd=2, col=c("navyblue"),
       horiz=F, bty="n", cex=.9, y.intersp=.8)

plot(1:8, basic.fragmentation[,1], type="l",
     ylim=c(0,.65), xlim=c(1,8), col="navyblue",lwd=3, lty=2,
     xlab="Tenure (years)", ylab="Prob. of political exit", 
     main="PANEL IV \n Simulated Effect of Expected Turnover \n when there is Fragmentation")
lines(1:8, basic.fragmentation[,2], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.fragmentation[,3], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.int[,1], col="navyblue",lwd=3)
lines(1:8, basic.int[,2], col="navyblue",lwd=.7, lty=2)
lines(1:8, basic.int[,3], col="navyblue",lwd=.7, lty=2)
legend("bottomleft", c("Expectation of Turnover", "No Expectation of Turnover"),
       lty=c(1,2), lwd=2, col=c("navyblue"),
       horiz=F, bty="n", cex=.9, y.intersp=.8)



###################################
# Appendix
###################################

# PARAMETRIC MODEL USING MAJORITY (1-fragmentation) AS MEASURE OF FRAGMENTATION
param.main.majority <- multinom(motive ~ tenure + t2 + t3 + majoritydesign + turnovert2 
                       + inter_fyturno, data = data)
param.main.complete.majority <- multinom(motive ~ tenure + t2 + t3 + majoritydesign + turnovert2 
                                + inter_fyturno + term_limit + poplog + head + lemas, data = data)


# PARAMETRIC WITH PROVINCE FIXED EFFECTS
#prov dummies
data$p1 <- ifelse(data$province==1,1,0)
data$p2 <- ifelse(data$province==2,1,0)
data$p3 <- ifelse(data$province==3,1,0)
data$p4 <- ifelse(data$province==4,1,0)
data$p5 <- ifelse(data$province==5,1,0)
data$p6 <- ifelse(data$province==6,1,0)
data$p7 <- ifelse(data$province==7,1,0)
data$p8 <- ifelse(data$province==8,1,0)
data$p9 <- ifelse(data$province==9,1,0)
data$p10 <- ifelse(data$province==10,1,0)
data$p11 <- ifelse(data$province==11,1,0)
data$p12 <- ifelse(data$province==12,1,0)
data$p13 <- ifelse(data$province==13,1,0)
data$p14 <- ifelse(data$province==14,1,0)
data$p15 <- ifelse(data$province==15,1,0)
data$p16 <- ifelse(data$province==16,1,0)
data$p17 <- ifelse(data$province==17,1,0)
data$p18 <- ifelse(data$province==18,1,0)
data$p19 <- ifelse(data$province==19,1,0)
data$p20 <- ifelse(data$province==20,1,0)
data$p21 <- ifelse(data$province==21,1,0)
data$p22 <- ifelse(data$province==22,1,0)
data$p23 <- ifelse(data$province==23,1,0)
data$p24 <- ifelse(data$province==24,1,0)


param.main.complete.provdummies <- multinom(motive ~ tenure + t2 + t3 + fragmentation + turnovert2 
                                            + frag.turnover + term_limit + poplog + head + lemas +
                                              p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 +
                                              p11 + p12 + p13 + p14 + p15 + p16 + p17 + p18 + p19 + p20 +
                                              p21 + p22 + p23 + p24 , data = data) # I exclude p1


# NO TIME CONTROLS
main <- multinom(motive ~ fragmentation + turnovert2 
                 + frag.turnover, data = data)
main.complete <- multinom(motive ~ fragmentation + turnovert2 
                          + frag.turnover + term_limit + poplog + 
                            head + lemas, data = data)
main.complete.provdummies <- multinom(motive ~ fragmentation + turnovert2 
                                      + frag.turnover + term_limit + poplog + head + lemas +
                                        p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 +
                                        p11 + p12 + p13 + p14 + p15 + p16 + p17 + p18 + p19 + p20 +
                                        p21 + p22 + p23 + p24 , data = data) # I exclude p1


## # NON PARAMETRIC MODEL
# Tenure dummies for nonparametric model
data$ten1 <- ifelse(data$tenure==1,1,0)
data$ten2 <- ifelse(data$tenure==2,1,0)
data$ten3 <- ifelse(data$tenure==3,1,0)
data$ten4 <- ifelse(data$tenure==4,1,0)
data$ten5 <- ifelse(data$tenure==5,1,0)
data$ten6 <- ifelse(data$tenure==6,1,0)
data$ten7 <- ifelse(data$tenure==7,1,0)
data$ten8 <- ifelse(data$tenure==8,1,0)
data$ten9 <- ifelse(data$tenure==9,1,0)
data$ten10 <- ifelse(data$tenure==10,1,0)
data$ten11 <- ifelse(data$tenure==11,1,0)
data$ten12 <- ifelse(data$tenure==12,1,0)
data$ten13 <- ifelse(data$tenure==13,1,0)
data$ten14 <- ifelse(data$tenure==14,1,0)
data$ten15 <- ifelse(data$tenure==15,1,0)
data$ten16 <- ifelse(data$tenure==16,1,0)
data$ten17 <- ifelse(data$tenure==17,1,0)
data$ten18 <- ifelse(data$tenure==18,1,0)
data$ten19 <- ifelse(data$tenure==19,1,0)
data$ten20 <- ifelse(data$tenure==20,1,0)
data$ten21 <- ifelse(data$tenure==21,1,0)
data$ten22 <- ifelse(data$tenure==22,1,0)
data$ten23 <- ifelse(data$tenure==23,1,0)
data$ten24 <- ifelse(data$tenure==24,1,0)
data$ten25 <- ifelse(data$tenure==25,1,0)

nonparam.main <- multinom(motive ~ fragmentation + turnovert2 + frag.turnover + ten2 + ten3 + 
                            ten4 + ten5 + ten6 + ten7 + ten8 + ten9 + ten10 + ten11 + ten12 +
                            ten13 + ten14 + ten15+ ten16 + ten17 + ten18 + ten19 + ten20 + ten21 + 
                            ten22 + ten23 + ten24 + ten25, data = data)

nonparam.main.complete <- multinom(motive ~ fragmentation + turnovert2 + frag.turnover + 
                                     term_limit + poplog + head + lemas + ten2 + ten3 + ten4 + 
                                     ten5 + ten6 + ten7 + ten8 + ten9 + ten10 + ten11 + ten12 +
                                     ten13 + ten14 + ten15+ ten16 + ten17 + ten18 + ten19 + ten20 + ten21 + 
                                     ten22 + ten23 + ten24 + ten25, data = data)

nonparam.main.complete.provdummies <- multinom(motive ~ fragmentation + turnovert2 
                                            + frag.turnover + term_limit + poplog + head + lemas +
                                              p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 +
                                              p11 + p12 + p13 + p14 + p15 + p16 + p17 + p18 + p19 + p20 +
                                              p21 + p22 + p23 + p24 + ten2 + ten3 + ten4 + 
                                              ten5 + ten6 + ten7 + ten8 + ten9 + ten10 + ten11 + ten12 +
                                             ten13 + ten14 + ten15+ ten16 + ten17 + ten18 + ten19 + ten20 + ten21 + 
                                              ten22 + ten23 + ten24 + ten25, data = data) # I exclude p1



# Tables for appendix

table.out <- function(model){
  m <- model
  return(xtable(cbind(summary(m)$coefficients[2,], 
                      summary(m)$standard.errors[2,])))
}


# Parametric with provincial dummies
table.out(param.main.complete.provdummies)

# No time controls
table.out(main) 
table.out(main.complete)
table.out(main.complete.provdummies)


table.out(nonparam.main)
table.out(nonparam.main.complete)
table.out(nonparam.main.complete.provdummies)






